Force distribution in a scalar model for non-cohesive granular material 
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We study a scalar lattice model for inter-grain forces in static, non-cohesive, granular materials, 
obtaining two primary results, (i) The applied stress as a function of overall strain shows a power law 
dependence with a nontrivial exponent, which moreover varies with system geometry, (ii) Probability 
distributions for forces on individual grains appear Gaussian at all stages of compression, showing 
no evidence of exponential tails. With regard to both results, we identify correlations responsible 
' for deviations from previously suggested theories. 

o' 

I. INTRODUCTION 

o 

Attempts to understand stress distribution in static, non-cohesive, granular materials have uncovered a rich structure 
t-J . of force chains that currently is poorly understood jjj . Although mechanical systems such as bead or granular packings 
subject to uniaxial compression seem straightforward, they harbor fundamental theoretical challenges concerning the 
connection between the micro and macro scales. Such systems are far from thermal equilibrium since thermal energy 
scales are negligible compared to the potential energies stored in elastic deformation of the grains. In addition, the 
force between grains is strongly nonlinear, as it vanishes identically when the grains are not in contact. 

In this paper, we report results and analyses of a toy model for the compaction of a granular material. Our model 
consists of "grains" placed on a square lattice and connected by vertical springs. The grains are constrained to move 
only in the vertical direction, so the model treats forces and stresses as scalar quantities. (See Fig. EL) Geometric 
disorder is introduced in the model by assigning to each spring an equilibrium length drawn randomly from a square 
distribution. Though we cannot hope to capture the true physics of tensorial stresses in this manner, we can study the 
£SJ ■ validity of arguments that have been applied to real granular systems and should equally well apply to our toy model. 
OS \ In particular, we compare the distribution of inter-grain forces generated by our model to the distribution predicted 
by the g-model of Coppersmith et al |Q , and we compare the stress-strain power law generated by our model to that 
expected on the basis of a mean- field argument. Our results provide a warning: in neither case are the theoretical 
predictions borne out by the numerics. 
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The structure of our model is shown in Fig. [I]. Grains are placed at the even-parity sites of a 2N X x (N y + 1) cubic 



lattice with lattice constant b. We refer to such a system as an "N x x N y " system. Excluding the sites on the lower 
boundary, it contains N x N y grains. Grains are constrained to move only in the vertical direction, as if sliding on fixed, 
frictionless wires (shown as thin lines in the figure). We let Uij denote the vertical displacement of the grain at lattice 
• • ■ site from its nominal height jb, with > for upward displacements. Between every two nearest-neighbor 
. ' grains there is a vertical spring. All of the springs respond linearly under compression, with spring constant k, but 
they produce no force under extension. The energy of the spring between grains i,j and i ± l,j — 1 is taken to be 

£± = ( 5 fc ~ s « ~ u i±W-i) if Ki - s tj - «i±i,j-i) < (!) 
tJ \ otherwise 

where sfj and s~- are independent, quenched random variables designating the different uncompressed lengths of 
the springs. We take sfj and s~- to be distributed uniformly over the interval (— e, e). Note that when s = 0, the 
uncompressed spring has length 6, the lattice constant. One may imagine the grains to have arbitrary shapes and 
the uncompressed spring length to be the vertical spacing between grain centers when the grains are just barely in 
contact. 

The system is taken to satisfy periodic boundary conditions in the horizontal direction. All grains at the top are 
constrained from above by a rigid "ceiling" and therefore satisfy < 0. On the bottom row, grains are constrained 
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by a rigid "floor" whose height U is an independent variable; thus u;.o > U. Let Uq denote the maximum value of 
U such that no springs are compressed. Typically Uq, which is a function of the random s^'s, is negative and of the 
order of N y e in absolute value. (It cannot be less than —N y e.) Below we measure the floor height in terms of the 
incremental variable A = U — Uq. 

Though we use the terms "floor", "vertical", etc., for convenience, there is no gravity in the model and there is 
complete statistical symmetry between the up and down directions. The algorithm discussed below for determining 
equilibrium configurations explicitly breaks this symmetry, though the configurations it produces do not. 



A. Meaningful parameters 

Our model apparently depends on the parameters e and fc, but these merely set the scales for distance and force. 
For the rest of this paper, numerical values for distances — spring-length variations s^-, grain displacements Uij, and 
boundary displacement A — will be quoted in units of e; and forces — individual forces /y and the total force F on 
the ceiling — in units of fee. 

Although the width of the distribution of spring lengths scales out of the problem, its shape could conceivably make 
a difference. Nevertheless, we do not expect the precise shape to be important as long as the probability of having 
very large uncompressed lengths decays sufficiently rapidly. In any event, we consider only the square distribution. 

Let us contrast the alternative of random spring constants with the random equilibrium lengths we study here. If 
the fc's were random but all equilibrium lengths were equal, then as A increased all the springs would begin to be 
compressed simultaneously at U = and the system would be perfectly linear for all A > 0; the force- vs.-displacement 
law would be F(A) = KA for A > and F(A) = for A < 0, where A' is a constant determined by the random fc's. 
By contrast, with random spring lengths we capture some of the geometric disorder of real granular materials, which 
inevitably introduces nonlinearities associated with the formation of new contacts during compression. 



B. General features of equilibrium configurations 

We will refer to a spring under compression as an "active bond" , and to a site connected to such a spring as an 
"active site." An equilibrium configuration for A positive, but not too large, has a mixture of active and inactive 
sites and bonds. For a given set of s^-'s, the equilibrium configuration of active sites and bonds is unique, although 
the exact positions of inactive sites are not determined. This uniqueness follows from the strict convexity of the total 
energy as a function of the lengths of active bonds. 

Fig. H shows three equilibrium configurations with the same set of s^-'s on a 40 x 40 lattice, for different amounts of 
compression. Only the active bonds are shown, with line thickness proportional to the force transmitted. Also shown 
are two bars representing the ceiling and the floor. 



3 shows a typical force curve for a 40 x 40 
ementary considerations reveal the limiting 



Let F(A) be the total force exerted by the system on the ceiling. Fig. 
system, together with the fractions of sites and bonds that are activated. E 
behaviors of F for large and small A. For very small A, the force on the ceiling is due to a single chain of activated 
springs that reaches from top to bottom. This chain, which follows the directed path through the lattice that 
yields the highest value of Ylj=i s y'i produces a force linear in A: 

F{A) = —A (A small). (2) 

Ny 

When A is very large, all of the springs in the system are compressed. The system is then again entirely linear: 
because all of the springs have the same stiffness, 

IN 

F(A) = ^-A + C (A large), (3) 

where C is a (negative) constant that depends on the random s^-'s and is difficult to compute. 

In general, as A is increased from zero, additional chains are activated and F(A) increases. In fact it can be shown 
that F must increase monotonically with A. (This is why C must be negative.) Surprisingly, however, the slope 
of -F(A) may decrease with increasing A. Occasionally bonds that are active at small values of A become inactive 
for larger A, thus reducing the effective spring constant of the network. The simplest example of this behavior is 
illustrated in Fig. ^|. The figure shows a sequence of configurations in which the central bond is originally active, but 
becomes inactive due to the effective stiffening of other parts of the network when additional bonds are activated. 
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III. THE NUMERICAL ALGORITHM 



The algorithm we employ for generating equilibrium configurations of our model relies heavily on the fact that for a 
given set of active bonds, the system response is linear. By solving for the times at which inactive bonds are activated 
during compression, the algorithm generates the entire curve F(A), starting from A = and explicitly visiting every 
configuration. To construct the initial, uncompressed configuration, we begin with Uifl — Uq for all i and march up 
layer by layer applying the rule 



max 



{(Ui-ij-i + s y -), (ui+ij-i + at)} ( 4 ) 



This configuration can be thought of as the packing that would result from an infinitesimal gravitational force acting 
on the system. (Strictly speaking, we set u^o = 0, march upward, determine Uq from the height reached, and translate 
all displacements accordingly.) 

A typical initial configuration is shown in Fig. |B[ Sites joined by line segments are separated precisely by the 
equilibrium length of the spring between them. Springs that are not long enough to connect their two sites are 
not drawn. Already one can see that there is a rich geometric structure hidden in this model, quite similar to that 
described by Roux et al || . The heavy lines in the figure highlight structures we call "trees" , consisting of all sites 
that are connected to the same site on the first layer. 

We refer to the distance between the upper end of an inactive spring and the site to which it will eventually connect 
as a "gap" . Note that sizes of the gaps reflect the structure of our initialization algorithm. We have moved all the 
grains to their lowest possible position, concentrating all of the small gaps that might be present in a chain into a 
single large gap at the top of a branch of a tree. 

To generate the force curve, we maintain a list of active sites and their active bonds, and we keep track of changes 
in the tree structure of the uncompressed chains as the floor is raised. For any fixed set of active bonds, the Ujj's 
increase linearly with A, since all of the active springs are linear. It is therefore a straightforward linear problem to 
solve for the rate, du/dA, at which each Uij advances. Advancing A through one stage involves solving for Uij for 
some arbitrary increase in A, for which we employ an iterative biconjugate gradient routine using the current 
configuration of itij 's as an initial guess for the solution. The solution is used to determine the rates of advance of all 
the active sites. The rate of advance of each inactive site is equal to that of the active site that supports it through 
an uncompressed chain. The rates are used to determine the value of A at which the network topology changes. The 
Ujj's are then updated according to the calculated rates, the list of active sites or the tree structure is updated, and 
the process is repeated. 

The updating of the network topology is necessitated either by the closure of a gap or by the breaking of a bond. 
The closure of a gap can generate one of two types of events. First, it can cause an additional set of sites to be 
activated, which increases the size of our active site list (and also increases the size of the linear problem that must 
be solved on the next iteration). Second, it can result in a "push" event, in which an inactive branch of one tree 
is simply transferred to another tree without becoming active. Pushes dominate the behavior of the system in the 
early stages of compression, becoming more and more rare as compression continues and the number of inactive sites 
decreases. 



The breaking of a bond, as mentioned in Section [IB. is the final possibility for changing the network of active 
bonds. Breaks are generally rare events. During full compression of a 40 x 40 system, which undergoes approximately 
2000 events (pushes, chain additions, breaks), there are typically about 5 breaks. 

The algorithm reaches completion when all bonds have been activated. As explained above, further compression 
would be homogeneous with all /^'s increasing at precisely the same rate. To compress one configuration on a 40x40 
lattice, generating one complete force curve, requires approximately one hour of computation on a typical 200MHz 
workstation. We have accumulated data for various lattices with at most a few thousand sites. Future efforts to 
optimize the algorithm should permit investigation of substantially larger lattices. 

Despite our reliance on linear methods for evolving the system between gap events, it must be emphasized that 
the absence of tensile forces introduces a strong nonlinearity for larger increments in A. An alternative approach 
to following the entire evolution of the system would be to solve directly for the configuration of the system at an 
arbitrarily chosen value of A by minimizing the nonlinear energy function of Eq. ([!]). Such an approach might speed 
up the calculation even with linear springs, and it would be essential if the springs were nonlinear (e.g. Hertzian) 
under compression. 
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IV. THE MACROSCOPIC FORCE AS A FUNCTION OF DISPLACEMENT 



A. General remarks and a "mean-field" prediction 

Fig H displays numerically computed curves of F vs. A for a few rectangular systems of various sizes. At a cursory 
level, the force curves appear to have the form F{A) ~ A", which correspond to lines of slope v on our log-log plots. 
Interestingly, this exponent depends strongly on system geometry. 

Closer inspection of Fig. ^| reveals both expected and unexpected behavior. For any particular realization of the 
disorder, there is a region, sometimes substantial, on the log-log plot for small values of A in which the force curve 
appears nearly linear because it is dominated by a single chain. There is also a crossover to linear behavior for large 
A since, as the fraction of active bonds approaches unity, the system approaches the linear limit described in previous 
sections. The intermediate regime is the one of interest, and the behavior there is rather complex. In many runs an 
exponent v > 1 in the intermediate regime can be readily identified. However, in other runs (see, e.g., Fig. |6|c), small- 
system statistics tend to obscure the phenomena. It appears that completely reliable measurements of the exponents 
in the various regimes will require bigger systems, beyond the reach of our current numerical codes. Nevertheless, we 
believe that the results for systems of a few thousand grains support the conclusions drawn below. (See Ref || for a 
discussion of a closely related model.) 

In studying force curves for bead packs, several authors have proposed a mean-field argument suggesting that for 
various granular systems one should observe v = a + 1, where a is the exponent of the single-contact force law. Thus 
for our model the mean-field theory predicts v = 2. (For the reader's convenience, we summarize a version of the 
theory in the Appendix. The treatment there is similar to that of References || and [||.) The argument is based on 
the assumption that the rates at which pairs of nearest neighbor sites approach each other may all be taken to be 
equal to the average rate of compression. 

In Subsection C below we present numerical results for our model showing that v can be significantly less than 2 for 
some system geometries and in Subsection D we seek to explain the failure of the mean-field argument. In our view, 
this failure casts doubt on the argument's applicability to real granular systems. Useful supplementary information 
is covered in Subsections B and E. 



B. Analytic results for limiting cases 



It is instructive to consider two cases for which the force curves can be explicitly calculated. 

Case 1: N y = 1. This system consists of only a single layer of random-length springs, as shown in Fig. ^a. Assuming 
the equilibrium spring lengths to be uniformly distributed in a finite interval, one immediately obtains F oc A 2 . In 
this case the mean-field argument (Appendix A) is exact, as all gaps do close at the same rate. 

Case 2: N x = 1. Because of the periodic boundary conditions, this case is equivalent to a single column of grains 
without periodic boundary conditions, as shown in Fig. ^|b. Again we take the distribution of equilibrium lengths to 
be a square distribution of width 2e. In this case it is more convenient to fix the force F and compute the displacement 
A(F) since the compressive stress on each grain must be the same. Let us decompose the displacement 

A = 5>„, (5) 

71=1 

where S n is the displacement of the nth spring. When F = 0, we have 6 n = for each n, the longer spring at each 
level being just at the threshold of compression. The growth of S n with F depends only on a n , the difference between 
the equilibrium lengths of the two springs in the nth layer. Specifically, 

d6n = f i if < F < ka n 
dF \ A. if kan < F. 

Note that the a„'s are independent and all have the probability density P(a) = (2e — a)/2e 2 for < a < 2e. 
The expected value of dA/dF is defined as 

Ny r 2e 

= £/ ^P(a n )da n . (7) 

n=l"' 
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Using the fact that all the terms in the summation are equal and performing the integration, we obtain (dA/dF), 
which in turn may be integrated to yield 

tX\ I f-lf 2 + llf 3 if < / < 2 

W= \|/ + | if 2 < /, < 8 > 

where S = A/N y e and / = F/ke are the nondimensionalized displacement and force per layer, respectively. Eq. (|^) 
shows that, in this case, there is no simple power-law behavior in the large N y limit. Note that Eq. (|^) is approximately 
linear near / = (i.e., limy^ + M( m ^)/^( m /)] = l)i an d it is linear near and beyond / = 2. If we attempt to identify 
a single power law for intermediate values of /, the natural choice is the derivative d(]xiS)/d(]xif), evaluated at the 
point where this quantity is most slowly varying; i.e., where d 2 (ln8)/d(lnf) 2 vanishes. This method yields / ~ 8 y 
with v = 1.27..., where we have returned to displacement as the independent variable. 

The analysis of Case 2 shows that it is possible for the force curve to exhibit a region corresponding to a power less 
than 2. It also shows that the emergence of a true power law should not be taken for granted in these systems. We 
find, however, that in sufficiently wide systems, a power law does arise (cf. discussion below). 



C. Numerical results 



Short, wide systems: Fig.|]a shows two typical force curves for an 80 x 20 system. It appears here that there is a 
regime in which v > 2, followed by the expected crossover to v = 1 at large A. This observation lends some support 
to the mean field argument and is consistent with the claim of Gilabert et al. [Q, who studied the electrical analogue 



of our model. (See Section IV E below.) 

Tall, narrow systems: In a 20 x 150 system, the exponent v appears to shift noticeably. The best fit to the curve 
shown in Fig. [j]b is v = 1.7 in the intermediate regime of interest. The error in this measurement is estimated to be 
±0.1 on the basis of fits made with different choices for which points to exclude from the intermediate regime. The 
data clearly rule out v = 2. 

Roughly square system: An estimate of v = 1.9 is obtained from data on 40 x 40 systems. Fig. ^c shows data from 
two 40 x 40 systems. The apparent lack of consistent behavior is presumably due to finite size effects. Nevertheless, 
using data from 25 runs, an estimate of v — 1.9 can be obtained as shown in Fig. ||d. To obtain the dotted line in 
the figure, the data in the intermediate regime were fit to a power-law; data from the initial linear regime and data 
from the two largest forces in the figure (where we expect a crossover to linear growth) were excluded in making this 
estimate for the exponent v. 

We suggest that the variation of exponent with aspect ratio can be traced to the tree structure in the system. 
Specifically, we conjecture that if most of the force is transmitted within a single tree, then a smaller exponent will 
be observed, whereas if the force is spread over many trees, a larger exponent. We interpret the following numerical 
experiment as support for this conjecture. Starting with a 60 x 60 grid, we removed all sites lying outside the 90° 
cone emanating from the center of the bottom boundary. We then measured F(A) for this reduced system. In this 
geometry all the active sites at every stage of compression are connected to the floor at the same point, so all the 
active bonds are contained within a single tree. As shown in Fig. ||e the exponent v is close to 1.5, indicating that 
the behavior within a single tree has a different character than in a system with many trees. 

To explore this issue further, we measured the rate at which trees expand as a function of height in the initial 
configuration. (Roux et al. have investigated a similar, but not identical, model ||.) While building the initial 
configuration, we keep track of the root of the tree associated to each site, thereby counting the number T(z) of trees 
that survive up to layer z. Since only initial configurations are involved — not their subsequent evolution — rather 
large systems can be simulated. For a system of width 50,000 with data collected for up to 10 7 layers, we obtained 
an excellent power law fit T(z) « Az^ 1 , with 7 = 0.66 and A — 6.5 x 10 4 , with only a slight deviation for very 
small heights (z < 10) and very large heights (corresponding to T(z) < 5). Note that A = 1.3N X . (Incidentally, since 
T(l) = N Xl it is clear that the power law fit must be inaccurate for z of order 1.) Generalizing from this system to 
one of arbitrary width, we estimate that in the initial configuration of an N x x N y system, approximately 1.3iV x iV~ 7 
trees will reach the top layer. 

Let us speculate on possible consequences of this estimate. The rates of advance of all sites within a single tree 
are correlated, since compression of the bottom bond of the tree affects the rates of all the branches above it. We 
therefore expect the sequence of gap closings to depend upon the initial tree structure. Assuming that the system 
obeys a simple anisotropic scaling law, we conjecture that in large systems the exponent v is a function of N X N~^ 
alone. Limited support for this conjecture comes from two additional runs on 20 x 80 and 40 x 225 systems. In 
both cases, N X N~ ' 66 s» 1.1, and in both cases we observe power-law behavior with v sw 1.8, as shown in Fig. ||f. A 
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more rigorous test would require runs on significantly larger systems, for which a more efficient code is needed. Note 
that the 40 x 225 system was already too large for us to run to completion; the curve in Fig. ^f ends well before the 
crossover to the linear regime occurs. Note also that in the 40 x 225 system, a clean power-law regime extends over 
at least two decades in A and three decades in F, giving us some confidence that a true power-law regime does exist 
in large systems. 

A consequence of our scaling conjecture is that increasing the system size while keeping the aspect ratio fixed should 
result in exponents that eventually approach that of a wide, short system. As we have shown that v = 2 in the oo x 1 
limit (Case 1, Section IVB), we expect to observe the mean- field result v = 2 in very large, approximately square 
systems as well. 

D. Why the mean field argument fails 

The key assumption in the mean-field argument is that nearest neighbor sites approach each other at the mean 
compre ssio n rate. In terms of the probability density for inter-site distances, the assumption is expressed quantitatively 
in Eq. ( |A5j) of the appendix, which asserts that the shape and width of the distribution of inter-site distances are 
independent of the compression A = N y 8. As discussed in the appendix, in our model the probability density P(x, S) 
is well defined only for x > 0; i.e., for bonds whose springs are compressed. However, even restricting our attention 
to the active bonds, we find that our data are inconsistent with the above assumption. Specifically, in Fig. ||a we 
show the probability densities for the (nonzero) forces /, where / = max{fca;,0}, at various stages of compression of 
a 40 x 40 system, and it is clear that the distribution broadens as S increases. More quantitatively, in Fig. |]b we plot 
the widths of the best fits to the data by Gaussian distributions restricted to {x > 0}, for several values of S. 

The broadening of P(x,5), which occurs during the early and intermediate regime of compression, is a simple 
consequence of force balance at branching points of chains of active bonds. Consider a site at which three active 
bonds meet, forming a "Y" (either right-side-up or upside-down). Because there are no tensile forces, the force in 
the unpaired branch of the Y is greater than either of those in the paired branches. Moreover, as the Y is further 
compressed, the rate at which the force in the unpaired branch increases must equal the sum of the rates of increase of 
the forces in the paired branches. Since the larger force evidently increases more rapidly than the two smaller forces 
(the forces almost always being increasing functions of A), the distribution broadens. 

During the late regime of compression, the addition of active bonds will convert Y's to X's, for which any correlation 
between force and bond compression rate must be more subtle. Indeed, when all of the bonds are activated, all bonds 
compress at exactly the same rate and the force distribution simply shifts uniformly to higher forces - the mean 
field assumption becomes exact. Thus we expect to observe a broadening of the force distribution during early and 
intermediate stages of compression, with a rate that decreases as the density of Y's becomes smaller. This is precisely 
the behavior displayed in Fig. j8| where the width of the positive x portion of the distribution is seen to grow roughly 
as a power of A less than unity for small A and level off at high compression. A quantitative calculation of the rate of 
broadening would require a detailed understanding of the statistics of branching in the active network and is beyond 
the scope of this work. 

Let us argue that the broadening of the force distribution generically leads to an exponent v smaller than the mean 
field value. Suppose, for definiteness, that P(x, 5) has a sharp leading edge at xi(5); i.e., that P(x, 5) = for x > Xx(8) 
and lim^^^- P{x, 5) is bounded away from zero. The existence of Y's in the active bond network leads to an advance 
of the leading edge xi(S) that is more rapid than the advance of 6 itself. Taking into account the formation of new 
contacts and branch points during compression, and assuming simple asymptotic behavior, we expect that Xi(S) ~ 5@ 
with j3 < 1. H Thus, to lowest order in S, we find F — J* P(x,S)kx dx ~ 6 2/3 , and therefore v < 2. Corrections 
to this exponent should become significant for 6 of order unity, which is also the order of the width of P(x). This 
expectation is consistent with Fig. pi where it can be seen that the crossover to the linear regime for large A occurs 
for A > Ny/5, which is the last half decade in the plots. 

These qualitative results do not depend on the linearity (under compression) of the springs in our model, nor do 
they depend on the dimensionality of the model. They should also apply to the vertical forces in a system with 
horizontal degrees of freedom, providing only that the creation of new branch points is sufficiently common. Thus 
the mean-field argument seems problematic for physical bead packs. The analysis shows that new contact formation 
can lead to exponents smaller than the mean field value, an important point in light of other mechanisms proposed 
to explain experimental observations of this exponent. 
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E. Relation to random resistor networks 



The analogy between the elastic properties of a network of linear springs and the electrical properties of the same 
network of resistors has been exploited in numerous studies []l0|| . In order to clarify the relation between our model 
and others, and particularly with random networks near the directed percolation threshold, let us consider the analogy 
in some detail. We will see that there are good reasons to be skeptical of the applicability of percolation model results 
to our model, but there is also an intriguing numerical coincidence. 

There is a formal identity between our equations for mechanical equilibrium 

ftj = ~ k (ui,j ~ u»±i,j - sfj) (9) 

and the electrical equations 

$ = ^(y itj -v i±ld -v±). (io) 

In the latter equation, I^j specifies the current between sites that are joined by a resistance R in series with a battery 
generating a potential difference vf^. Note that the resistors' conductance plays the role of the springs' stiffness. The 
property that the springs function only under compression can be modeled in the electrical system by the insertion of 
perfect diodes, all directed "downward", in series with each resistor. Roux et al. studied this very model but did 
not report results for varying aspect ratios or single trees and did not study its relation to the g-model, which had 
not yet been introduced. 

A system with batteries is substantially more complicated than a simple, randomly diluted resistor-diode network. 
To relate the randomly diluted network at the directed percolation threshold to our model just beyond the initial 
linear regime, one must first assume that the network of active springs in our model has the same structure as the 
current carrying paths in a directed network of resistors placed at random on the lattice. One must also assume a 
relation between A of our model and the probability p that a bond exists in the percolation system. The natural 
assumption here would be that the probability p—p c is proportional to A, where p c is the critical value for percolation, 
since A = is the point where a single force chain (or current carrying path) first forms and on average the bonds in 
our model are compressed by an amount proportional to A []l0|| . In our system, however, the batteries play two roles. 
First, they generate potentials that affect the current distribution even when all the diodes are forward biased, which 
would correspond to a trivial, homogeneous state of the simple resistor network. Perhaps more importantly, however, 
they determine which diodes will be forward biased for a given applied potential difference across the whole network. 
This dynamical process of selecting the current carrying paths may or may not yield structures well-modeled by the 
random addition of resistive links in a percolation model. 

In spite of the difficulties in establishing a connection, it is interesting to compare our results for the power law 
obeyed by the stiffness dF/dA to the conductivity exponent obtained from the theory of directed percolation in 
randomly diluted resistor networks [ pd| , ^2[ . The conductivity exponent for directed percolation has been calculated 
both numerically and using renormalization group methods and appears to be approximately 0.7 ±0.1, |Tl] , |l^ ] which 
corresponds to a value of 1.7±0.1 for the exponent v. The fact that this agrees with our measurements on tall, narrow 
systems deserves further study. Other authors have investigated additional details of the statistics near threshold in 
our model and argued that the system is closely related to the percolation one. |13] 



V. STATISTICS OF FORCES ON INDIVIDUAL GRAINS 



We now consider the statistics of forces transmitted by individual springs, exploring, in particular, the relation of 
our results to the g-model of Coppersmith et al. [^j . Our model is simpler than the granular packings the q- model was 
intended to describe. However, our model would appear to be a better candidate for description by the g-model than 
the original granular packings, since we have removed the tensorial stresses from our system but still have contacts 
whose formation is governed by quenched randomness. 

In the q- model in two dimensions, qij refers to the fraction of the force on site from above that is transmitted 
to its neighbor at (i — 1, j — 1), the complementary fraction 1 — qi j being transfered to site + — One assumes 
that each qi j is a random variable that (1) is independent of the other g's and (2) has the same distribution 77(g) at 
every site, independent of the force supported by that site. 

Analytical studies of the g-model show that the distribution of forces supported by a single grain has an exponential 
tail at large forces whenever rj{q) is non- vanishing at q = 1. [p|fl4|] By contrast, the force distributions in our model, 
plotted in Fig. 0, show no evidence of exponential tails. At no stage during the compression does it appear that the 
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g-model distributions are a good match for the distributions we observe. In particular, consider the fourth curve from 
the left in Fig. ^, which is made at a compression for which almost all sites are active but there remains an appreciable 
fraction of inactive bonds. The distribution of (ftj's directly measured from our data in this regime (see Fig. ^|a) is 
reasonably represented by 

rj(q) = 0.1 [%) + (5(1 - q)} + 0.8 - 0.2 cos(27rg), (11) 

for which 20% of the g's are either or 1. The g-model would predict an exponential tail in the single grain force 
distribution. Q Even for the small systems we study here, the exponential tail would be clearly distinguishable from 
the rapid decay we observe. Fig. ^|b compares numerical results from our model and from a simulation of the g-model 
with 77(g) given by Eq. (pi]). For the q- model simulation, forces on the top row were chosen randomly from a uniform 
distribution on the interval (0,2). 

To understand the discrepancy between our results and the predictions of the g-model, we examine our data as 
regards assumption (2) above: i.e., that the distribution of q at each site is independent of the force supported by 
that site. Fig. |9jc shows the distribution of q's obtained from equilibrium configurations of our model corresponding 
to the same conditions as in Fig. ^|a, but separated according to the force supported by the site. Different symbols 
in the plot indicate different levels of force as described in the figure caption. It is clear that the contribution to 
r)(q) from larger forces is peaked more strongly about 1/2 and has little weight near or 1, which explains why 
the g-model prediction fails for this system. As a quantitative measure of the correlation, we have computed the 
covariance C qw = ((qtj — \) 2 {wij — (wtj))) for our model, where Wij = + and the averages are performed over 
space for one realization of a 40 x 40 system. (The square in this definition is necessary because left-right symmetry 
guarantees that a correlation function linear in q — 1/2 would vanish.) C qw vanishes identically in the g-model since 
qi.j and Wij are independent in that model. As shown in Fig. |^d, C qw develops a significant negative value when the 
system is under compression, indicating that larger forces are associated with g's closer to 1/2. The heavy dot in the 
figure indicates the point corresponding to the data in parts (b) and (c). 

VI. CONCLUSIONS 

We briefly summarize the results of our study of the scalar model, discuss its generalization to three dimensions, 
and finally draw two conclusions concerning the implications for real systems or more realistic models. 

In the context of a toy model, we have tested arguments that have been applied to stresses in static, non-cohesive 
granular materials. Further study of the model is needed, especially simulations of larger systems, but already two 
important facts have been established. (1) Correlations in the stress configuration are responsible for substantial 
effects, both at the level of single-grain forces and that of macroscopic stresses. (2) Depending on the (properly 
scaled) aspect ratio, new contact formation may play a decisive role in determining the macroscopic stress-strain 
relationship, with sufficiently tall systems showing power-law behavior with a nontrivial exponent. Regarding (1), 
we have identified two important effects: (i) the distribution of forces broadens under compression because of force 
balance constraints at sites where three active bonds meet; and (ii) larger forces tend to divide more evenly between 
supporting grains as a result of the dynamical process that determines the structure of the active bond network. 
Theories that neglect these correlations fail when applied to our model. 

One may wonder whether qualitatively new features might appear in a 3D generalization of our model in which the 
vertical direction is taken to be the 111 direction of a simple cubic lattice. As a preliminary check, we have measured 
the statistics of trees in the initial configuration and observed their general morphology. The number of surviving 
trees decays as z -12 and the trees remain relatively compact. In other words, the diameter of surviving trees grows as 
z o.60 ( com p are d to z - 66 in two dimensions) and the branches of separate trees do not become heavily entangled with 
each other. On this basis, we conjecture that the essential physics of the stress-strain relation will not be qualitatively 
different in three dimensions. In particular, we expect a variation in the exponent v with the scaled aspect ratio. 

We draw two general conclusions. First, it appears that the mechanism for the broadening of the vertical force 
distribution should be present in the full tensorial problem as well, and therefore should be re-examined as a possible 
explanation for the observation of smaller exponents than those derived from the mean-field argument. Experiments 
measuring the dependence of the exponent v on aspect ratio would be especially interesting. Second, the success of 
the g-model in predicting the distribution of single-grain forces in bead packs stands in need of explanation. The 
g-model has already been criticized for its neglect of tensorial stresses |lq| . One may have expected it to provide 
a better description of the stresses in a scalar model of the type we have studied, but this is not borne out by our 
results. 

We thank Sue Coppersmith and Long Nguyen for helpful conversations, and also Scott Zoldi for his suggestions. 
This work was supported by NSF through grants DMS-98-9803305 (DGS and MGS) and DMR-94-12416 (JESS). 
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APPENDIX A: THE MEAN-FIELD ARGUMENT PREDICTING v = 2 



In this appendix, we summarize a mean-field argument that was developed to describe forces during the compression 
of a bead pack between parallel plates (See, e.g., Refs. 

We define the overlap between two adjacent beads as (do — d) where d is the distance between their centers and do 
is the nominal distance at which the beads touch but exert no force. Note that because of the minus sign, the overlap 
is positive for beads in contact, and negative for beads not in contact. 

We introduce the random variable, X to be the overlap between two randomly chosen adjacent beads. Thus the 
sample space for X includes both the choice of configuration (spring lengths in our model), and the choice of a pair 
of adjacent beads. Of course X also depends on the displacement A of the floor plate. (As in Section I, we normalize 
so that nonzero forces start at A — 0). We rescale the independent variable, defining 8 = A/N y where the bead pack 
is N y layers thick. Let P(x,S) be the probability distribution for X. The average force per spring, a, is given by; 



a(6) 



f(x)P(x,S) dx, 



where f(x) = max{fca:,0}. Since 

f(x) = for x < 0, 
the integral may be restricted to (0, oo). Of course, by our normalizations, 

P(x,0) = for x>0, 



(Al) 



(A2) 



(A3) 



so tr(0) = 0. The rescaled variable 5 = A/iV a equals the average change in the overlap X resulting from motion of 
the floor; thus, 



xP(x, S) dx = 5 + So, 



(A4) 



where So is the average overlap when S — 0. 

The first main assumption of the mean-field argument is a considerable strengthening of ( A4): One assumes that 
the shape and width of the probability distribution are independent of S, i.e., 



P{x,S) =P(x-S,0). 

Substituting (AE) into ( [Al] ) and rescaling (A^,A3|), we deduce that 



cr = k [ (x + S)P(x, 0)dx, 
J-s 



which is equivalent to Eq. (1.2) of Ref. ||. 

The second main assumption of the mean-field argument is that P is continuous in both arguments and 

lim P(x,0) = C, 

x— >0- 

where C is a finite positive constant. Combining the two assumptions, we find that 



(A5) 



(A6) 



(A7) 



k (x + S)[C + 0(x)} dx, 



-5 



(A8) 



which immediately yields v = 2 . Under assumption ( A5), the only way to accommodate a value of v less than 2 is to 
posit that the limit in Eq. (A7) is infinite. 

Note that in our model, the distribution for X is not uniquely defined for X < sin ce in active grains are free to 
relocate slightly. Despite this ambiguity, in section IVD we are able to test hypothesis (A5) by focusing on positive 
values of X. 
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Figures 




FIG. 1. The scalar model. Diamonds represent rigid grains that slide on frictionless, vertical wires. Squiggly lines represent 
springs with identical spring constants but differing equilibrium lengths. Grains cannot penetrate the ceiling or floor. 
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FIG. 2. Configurations of active bonds in a 40 x 40 system at three different stages of compression. 
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FIG. 3. A typical force curve (solid line) for a 40 x 40 system together with curves showing the fraction of active sites (dotted) 
and the fraction of active bonds (dashed). 
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FIG. 4. A sequence of configurations in which the central bond is originally active, but becomes inactive after other bonds 
become active. Consequently the total stiffness of the network decreases as the system is compressed. Thick (thin) solid lines 
represent active (inactive) springs with s = 0. Thick (thin) dashed lines represent active (inactive) springs with equilibrium 
length s = — e. Springs not drawn are assumed to have small enough equilibrium lengths that they are never active in this 
sequence. dF/dA is smaller for the third configuration than for the second, although A itself has increased. 




FIG. 5. A typical initial configuration. All springs shown are at exactly zero compression. Springs with equilibrium lengths 
too short to connect their pair of sites are not drawn at all. The heavier lines are guides to the eye, highlighting two of the 
trees in this configuration. 
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FIG. 6. Plots of F(A) for several different system sizes. Line segments with labels indicating their slopes are added as 
guides to the eye. Labels in the lower right corners indicate system size, (a) A short, wide system showing an exponent near 
2. (b) A tall, narrow system showing an exponent near 1.7. (c) A square system showing apparently different behavior for 
different realizations of the quenched random equilibrium lengths, (d) Data from 25 runs of a 40 x 40 system. For clarity, 
points are plotted only at 10 discrete values of F. (Note the different horizontal scale from (c)). The dotted line through the 
data is a power-law fit to the points shown, excluding the two highest forces. The exponent is 1.9. (e) A single tree showing 
an exponent of 1.5. (f) Two different system sizes with equal values of N x N~ ose ', both showing an exponent of 1.8. The open 
circles represent one completed run of a 20 x 80 system. The filled circles represent one partially completed run of a 40 x 225 
system. 
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(a) 




FIG. 7. Limiting models for which F(A) is easily computed analytically: (a) an infinite row of springs with random lengths 
(N v = 1); (b) an infinite column of single grains with two springs in every row (N x = 1). 
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FIG. 8. Distribution of spacings between nearest neighbor grains for several stages of compression of a 40 x 40 system, 
(a) From narrowest to broadest, the curves represent active bond densities of 0.2, 0.4, 0.6, 0.8 and 1.0, which correspond to 
ensemble averaged total forces of 0.85, 4.3, 11.2, 27.9, and 116.1 xfce. The data represents averages of 25 runs, (b) The width 
of the fitted Gaussians as a function of A on a log-log scale. Horizontal error bars indicate the range of values of A obtained 
for a given density of active bonds in different runs. The thin lines in the left plot are the Gaussians used for the maximal and 
minimal width estimates at each stage. The line of slope 1.0 on the right is a guide to the eye. 
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FIG. 9. Distribution of q values and comparison with predictions of the q-model. (a) The frequency distribution of q's 
averaged over 25 configurations of a 40 x 40 system. The different symbols indicate different stages of compression, with larger 
values at q = 1/2 corresponding to larger compressions. All distributions are normalized to unity, but the points at q = and 
q — 1 are off scale in all cases except for the fully compressed one. (b) Comparison of the distributions of total force supported 
by a single grain in our model and in the q- model with a similar frequency distribution of q values. (See text for details.) 
The open circles are averages over the same 25 configurations used to generate the open circles in (a). The heavy dots are 
averages over 100 configurations of the q-model. (c) Separation of the q distribution from one curve in (a) into components 
corresponding to different levels of supported force. The data in (c) correspond to the open circles in (a). The bins used are 
(0, W) (open circles), (W, 2W) (open triangles), (2W, 3W) (filled triangles), and (3W, oo) (filled circles), where W is the average 
force supported by a single grain, (d) The correlation function C qw = ((qij — |) 2 (wij — ( w ij))}- The data shown are from one 
realization of a 40 x 40 system. The heavy dot indicates the point corresponding to the data in (b) and (c). 
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